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Abstract 

A one dimensional fractional diffusion model with the Riemann-Liouville fractional 
derivative is studied. First, a second order discretization for this derivative is pre- 
sented and then an unconditionally stable weighted average finite difference method 
is derived. The stability of this scheme is established by von Neumann analysis. Some 
numerical results are shown, which demonstrate the efficiency and convergence of 
the method. Additionally, some physical properties of this fractional diffusion sys- 
tem are simulated, which further confirm the effectiveness of our method. 
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1 Introduction 



Recently, a large number of applied problems have been formulated on frac- 
tional differential equations and consequently considerable attention has been 
given to the solutions of those equations. Fractional space derivatives are used 
to model anomalous diffusion or dispersion, a phenomenon observed in many 
problems. There are some diffusion processes for which the Pick's second law 
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fails to describe the related transport behavior. This phenomenon is called 
anomalous diffusion, which is characterized by the nonlinear growth of the 
mean square displacement, of a diffusion particle over time. The anomalous 
diffusions differ according to the values of a, where a is the order of the frac- 
tional derivative. Some works providing an introduction to fractional calculus 
related to diffusion problems are, for instance, [2,6,11,12,28,29]. In this work 
we will be interested in the anomalous diffusion, called supperdiffusion, for 
1 < a < 2 and experimental evidence of this type of diffusion is already 
reported in several works [1,7,13,14]. 

Fractional derivatives are non-local opposed to the local behaviour of inte- 
ger derivatives. Therefore, different challenges appear when we try to derive 
numerical methods for this type of equations. Numerical approaches to dif- 
ferent types of fractional diffusion models are increasingly appearing in lit- 
erature. We can found recent work on numerical solutions for the fractional 
diffusion equation describing superdiffusion [5,9,10,19,15,27,16] and also for 
several transport equations including this type of diffusion [18,25,31]. Some 
other works consider subdiffusion, which is represented by a time fractional 
derivative of positive order and less than one [3,30]. However, the challenges 
for these equations are different from the ones that arise when we consider a 
space fractional derivative of order 1 < a < 2. 

Numerical methods, for models with superdiffusion, have been obtained with 
mathematical techniques which do not necessarily consider a second order 
discretization for the fractional derivative to achieve second order accuracy. 
In this work, we present a second order approximation for the fractional Rie- 
mann-Liouville derivative of order a, 1 < a < 2. This approach uses some of 
the tools described in [4,8] and also applied in [26] to derive an approximation 
for the Caputo fractional derivative defined in bounded domains. Here, we 
consider the Riemann-Liouville fractional derivative in an unbounded domain 
and its discretization is represented by a series instead of a finite sum. We 
prove the order of consistency of this discretization is second order. 

A weighted average finite difference r-scheme is considered, for r G [1/2,1], 
which includes the Crank-Nicolson method (r = 1/2) and the back forward 
Euler method (r = 1). The consistency and stability of the r-scheme are 
established and we prove the r-scheme is unconditionally stable. Also for r = 
1/2 we have second order accuracy in time and space as expected. 

Consider the one-dimensional fractional diffusion equation [1,7,16] 



on the domain a; G IR, where 1 < a < 2 and d{x) > 0, subject to the initial 



du 
dt 



(x, t) = d{x) 



dx' 
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condition 



u{x,0) = f{x), X eJR 



(2) 



and to the boundary condition 



u(x,t) = as |x| — )■ cxD. 



(3) 



The usual way of representing the fractional derivatives is by the Riemann- 
Liouville formula. The Riemann-Liouville fractional derivative of order a, for 
X G [a, b], — oo < a < 6 < oo, is defined by 



where r(-) is the Gamma function and n = [a] + 1, with [a] denoting the 
integer part of a. 

The function u{x,t) under consideration, that is, which is solution of (1), 
should be such that the corresponding integral (4) converges. If the function 
u{x, t) vanishes at infinity, as assumed when we impose the boundary condition 
(3), we have absolute convergence of such integrals for a wide class of functions 
[24]. However, these functions do not necessarily need to vanish at infinity 
and we can found under which conditions these integrals converge in [24] 
(section 14.3). There are very complete works about the fractional calculus 
[17,20,21,22,24], where the theoretical properties of this type of derivative are 
studied in detail. 

Another way to represent the fractional derivative is by the Griinwald-Letnikov 
formula, that is. 



The Griinwald-Letnikov approximation is often used to numerically approxi- 
mate the Riemann-Liouville derivative and it was the first algorithm to appear 
for approximating fractional derivatives [21,22]. However, this approximation 
has consistency of order one and also very frequently numerical approximations 
based in this formula originate unstable numerical methods and henceforth a 
shifted Griinwald-Letnikov formula is used [16,18]. 

The plan of the paper is as follows. In section 2 we derive a numerical ap- 
proximation for the Riemann-Liouville derivative. The full discretisation of 




(n - 1 < q; < n) (4) 
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the fractional diffusion equation is given in section 3, where a weighted finite 
difference method in time is apphed with the weight r G [1/2, 1]. In section 
4 we prove the convergence of the numerical method by showing consistency 
and stability. In the fifth section we present numerical results which confirm 
the theoretical results and in the last section we give some conclusions. 



2 The numerical method 



In this section we present a numerical approximation for the Riemann-Liouville 
derivative and also the numerical method that gives an approximate solution 
to the fractional diffusion equation. 



2. 1 Approximation of the Riemann-Liouville derivative 



Let us consider the Riemann-Liouville derivative [21,22], that is, 
d'^u , , Id 



We define the mesh points Xj = jAx, j & 1^ where Ax denotes the uniform 
space step. For a fixed time t, let us denote 



laix) 



J ui^,t)ix-0'~''d^- (7) 



First, we do the following approximation at xj 
^2 1 

—X^{Xj) - [^a{Xj-l) - 2Xa{Xj) + Xa{Xj+l)] . 



For each Xj we need to calculate Xa{xj). 

We compute these integrals by approximating u{^,t), at a fixed instant t, by 
a linear spline Sj{C,), whose nodes and knots are chosen at Xk, k = . . . , j — 1, j , 
that is, an approximation to Xa{xj) becomes Ia{xj) defined by 



■^3 
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The spline Sj{^) interpolates the points {{xk,t) : k < j} and is of the form 
[23] 



fc=— oo 



(9) 



with Sj^kiOy sach interval [xk-i,Xk+i], for A; < j — 1, given by 



1 



Xk+1 - ^ 
Xk+1 ~ Xk 



, Xk^l <^<Xk 



Xk<i< Xk^i 



otherwise, 



(10) 



and for k = j, 



Xj 



otherwise. 



(11) 



From (8) and (9), 



k=—oo 



Xk+i 



(12) 



We have that 



^k+i 



^ - Xk-l 



Ax 



^k+i 



Ax 



Xk-l 



Xk 



Ax 



2-a 



(2 - a)(3 - a) 



(13) 
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where the a, fc are such that, 



(j -k + 1)3-" - 2(j - + (j -k - 1)3-", k<j-l 



1, 



(14) 



k = j. 



Therefore, 



la (Xj 



(2-a)(3-a) 



(15) 



and an approximation for — — ■Xq,(x,), is given by, 

OX'' 

1 



Ax2 



[/c,(xj_i) - 2Ia{xj) + Ja(xj+i)] 



(16) 



that is. 

Ax-" 



(2 - a){3 - a) 



n(xfc,t)aj_i,fc - 2 ^ n(xfc,t)aj_fc+ ^ u{xk,t)aj+i^k 

k=—oo k=~oo k=—oo 



Let us assume there are approximations U" := {U'j'} to the values u{xj,tn) 
where tn = nAt, n > and At is the uniform time-step. 



We define the fractional operator as 
1 



S 



i+1 



(17) 



where 



cij,k = «i-i,fc ~ '2aj,k + o-j+i,k, k < j — 1 
Qjj+i = aj+ij+i. 

Therefore, an approximation of (6), for t = tn, can be given by 



(1^ 



Ax"' ' 
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We can also write the fractional operator (17) as 

1 oo 

^^U^ = Y^^E^^^.-mUf-m- (19) 



Remark: Note that for a = 1 and a = 2 the coefficients (18) are such that 
Qj^k = 0, for k < j — 1. For a = 1, qj^-i = —1, Qjj = 0, qj^+i = 1 and for 
" = 2, qj,j-i = 1, Qjj = -2, Qjj+i = 1. 

Remark: The series (19) converges absolutely for each 1 < a < 2 and for 
every bounded function u{x,t), for a fixed t. This result is a straightforward 
consequence of some results given in section 3 about the convergence of the 
series of the qjj-m- 

In this section we have considered a linear spline to approximate the inte- 
gral representation of the Riemann-Liouville derivative with the purpose of 
obtaining a second order approximation. In the next section we describe the 
full discretisation of the differential equation. 



2.2 Weighted average finite difference methods 



We discretize the spatial a-order derivative following the steps of the previous 
section. The discretization in time consists of the weighted average discretiza- 
tion. 

We consider the time discretization < t„ < T. Additionally, let dj = d{xj), 
pu _ p(^Xj,tn)- For the uniform space step Ax and time step At, let 

a _ djAt 



From equation (1) we can arrive at the explicit Euler and implicit Euler nu- 
merical methods, respectively 

^ ^-^W+P^ (20) 



At Ax"- 



jjn+l _jjn , 

^ = - -^s^m^' + p-+\ (21) 



At Ax"- 



Let (20) multiplies (1 — r) and (21) multiplies r. We obtain the following 
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weighted r-scheme 

^n+i _ jjn ^ _ ^)^^^n ^ ^^^f/^+ij + + (1 - r)Atp^", (22) 



where r G [1/2, 1]. 

Note that for a = 2, the operator (17) is the central second order operator 
that is, 

5at/; = f/;Vi-2f/; + f/7-i. 

We have the following numerical method 

(l - = (l + (1 - r)/.J5,) f/; + Atpf ^ (23) 

where 

= rp"] + - r)'p']^\ 



3 Convergence of the numerical scheme 



In this section we prove the convergence of the numerical method by showing it 
is consistent and von Neumann stable. First, we start to study the consistency 
of the numerical method and lastly we present the stability results. 



3. 1 Consistency 



In the beginning of this section, for the sake of clarity, we omit the variable t 
and we denote the partial derivative of m in x of order r by u^'^\ 

Lemma 1 Let u e C(^)(1R). For ^ G 

1 3 ^ 

W(0 -Sj,fc(0 = f IIw^''n04,r(0 - TyM^^n%)4,r(0> % e [Xk-l,Xk], 

^- r=2 ^■ 

where 



|4,,(0| < Ax^ 

Proof: For ^ G 

u{0 - SjAO = u{0 Aa^^*^^'^'"^'* Ax~^*^^^'*' 

Using Taylor expansions, we obtain 
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1 1 



r=2 ^■ 

where lk,r{0 functions which depend on Ax and Xk, given by 



r 



{xk-0^ + J2 ' \{xk-0'^\-'^y-'^x''-'-'- (25) 



r=0 



It is easy to conclude that |/fc,r-(OI — '^x^, for ^ G [xfe„i,Xfc]. □ 

Theorem 2 Order of accuracy of the approximation for the fractional deriva- 
tive): Let u G C^^\]R) and such that u^'^\x) = 0, for x < a, being a a real 
constant. We have that 



|ei(x,)| < CiAx^ \e2{xj)\ < CsAx^, 
anc? Ci and C2 are independent of Ax. 

Proof: It is straightforward to prove that we have 

= r(2 - a) Ax^ [^"(^J-i) ~ '^Mxj) + Mxj+i)] + 
where ei(xj) = C(Ax^). 
Let us define the error Es{xj), such that, 

Ja(Xj„i) - 2Xq(Xj) +Xa{Xj+i) = Ia{Xj-i) - 2Ia{Xj) + /^(Xj+i) + Es{Xj). 

We have 



d°'u 1 1 

d^^""'^ ^ r(2 - a) Ax^ [^"(^^-1) - + ^"(^^+1)] 



that is 
where 

We are now going to compute the error Es{xj). We have 



fc = -00 

i+1 



Taking in consideration the previous lemma, let us denote 
^ 1 

r=2 ' ■ 



where Ej.{ defined as follows. For r = 2 and r = 3, 



i+1 7= 



i+i 
+ ] 

k — 00 Xk-l 

and for r = 4 



j 

k = -00 
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For r = 2, 3 by changing variables, we obtain 



j 

k=-oo 
j 



A:=— oo 



^k-l 



that is, 



1 

k=—oo 



Let Xq = Na^x such that ^'•^^(x) = 0, for x < Xa- For r = 2 we have 



Xk 



E2{x,) = E / ikAO h^'^He + Ax) - 2««(o + u^'H^ - Ax)] (x, - o'-'^d^ 



k = -00 



Ax^ i 

— — E ^i^^H6)Ci,fc,2, 6e[Xfc-l,Xfc] 



k=Na+l 



where 



Since, by Lemma 1, 



xk 

xk-1 



Xk 

|c„fe,2| < Ax^ j [x.-if-'^di 

xk-1 



and 



we have 



^3 

hx,-if-^di=- (x,-x„) 

J Z — a 



2~a 



A^4 

\E2{xj)\ < ^ r||nW|U(x,- -x,)2-° 



2(2 -a) 



(29) 
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For r = 3, 



EsiXj) = J2 ^^iu^ i^ki) - M^^n^fc2))Cj, fc,3, 61,62 e [Xk-l,Xk\ 

k = Na + l 



and 



We have 



O A^4 

mx,)\ < r||«(^)||oo(x, -X.)^-" (30) 

(2 -«) 



Finally for r = 4, we bound each integral of (28) separately. For the first 
integral we have 



<Ax^||w(^)|U E / (^.-i-0'-"rfe 



Ax^ (4) 2 



2- a 

Therefore, since (a + hy < \a\^ + \h\^ for < p < 1, we have 
E / /M(0(^.-i-0'""^^e<^^ll^^'^l|oo((xi-a;.)2""+Ax2-") 



E ^^'H^fc) / M0(^.-0'-"rfe<;7^ll^^'^l|oo(a:,-xJ 



2-a 



Similarly, for the second integral we have 

Ax 

/M(^)(^i-^)^-"d^< ■ 

k=Na + l Xk-1 

and for the third integral 

E ^^'^(%) / /M(0(a^.+i-0'-"^^e<^^ll«^'^l|oo((x,-x„)2-"+Ax2-"^ 
Finally, we have 

OA 4 OA 6-a 

|£^4(x,)|<i^lk')||oo(x,-xj2-° + ^^^|M^)|U. (31) 
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Prom (29), (30) and (31) it is easy to conclude that the error Es{xj) defined 
by (45) is of order O(Ax^) and therefore the e2{xj) is of order 0{Ax^). 



□ 

Theorem 3 The truncation error of the weighted numerical method (23) is 
of order C(Ax^) + C(At'"^), where rrir = I, for r G (1/2, 1] and rrir = 2, for 
T = 1/2. 

Proof: Let u = u{x, t) be a solution to the fractional partial differential equa- 
tion and satisfying the conditions of the previous theorem. Note that the 
truncation error for the numerical method (23) is given by 

We have that 

- _ Q^^M ^ At dM^j,Q ^ ^^^^2)^ (32) 



At dt 2 dt 

and using the previous theorem we have 



.n_ du{xj^Q , At d'^u{Xj,tn) ^ ^l^^2^ fj d'^u{xj,tn+i] 
^ ~ dt +2 dt^ + ^ ^ { ' dx- 



_^,_,^d/-y^ + OiAx^))-pr 



Therefore 



_ du{Xj,tn) Atd'^u{Xj,tn) , du{Xj,tn) du{Xj,tn 

' ~ dt 2 dt^ ^ ^' dt ^ dt 

+0{At') + 0{Ax^) 



Finally, 



T^ = {\- r) At^^^^||l^ + 0{Ae) + O(Ax^) 



□ 
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3.2 Fourier decomposition of the error 



In order to derive stability conditions for the finite difference schemes, we 
apply the von Neumann analysis or Fourier analysis. Fourier analysis assumes 
that we have a solution defined in the whole real line. It is also applied to 
problems defined in finite domains with periodic boundary conditions since 
the solution is seen as a periodic function in IR. 

If is the exact solution let 

^1 = ^1 - ^" (33) 



be the error at time level n in mesh point j. To apply the von Neumann 
analysis we also consider dj locally constant, and we denote /i^ by /i". 

Considering the scheme (23) and inserting equation (33) into that equation 
leads to 

(1 - r^"5„) = (1 + (1 - r)/i"5,) El (34) 

The von Neumann analysis assumes that any finite mesh function, such as, 
the error will be decomposed into a Fourier series as 

N 

E]= f^y^'^'^^K j = -N,...,N, 

p=-~N 

where is the amplitude of the p-th harmonic and C,p = pir/NAx. The 
product ^pAx is often called the phase angle 9 = ^pAx and covers the domain 
[— 7r,7r] in steps of ir/N. 

Considering a single mode n'^e^^^, its time evolution is determined by the 
same numerical scheme as the error E^ Hence inserting a representation of 
this form into a numerical scheme we obtain stability conditions. The stability 
conditions will be satisfied if the amplitude factor k, does not grow in time, 
that is, if we have \n{6)\ < 1, for all 6. 

As we have seen the fractional operator can be written as 

-1 oo 

where the qjj-m are defined by (18). 

First we plot, in Figures 1-2, the coefficients qjj-m and then we give the 
properties that allow us to conclude this is a well-defined operator. 
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-1^ 
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a 



1.8 




Fig. 1. Coefficients (18): (a) qjj (b) Qjj-i 

0.014 




Fig. 2. Coefficients (18): (a) qj 



,j-m, 



m 



2, 3, 4; (b) qj,j-m, m = 5, 6, 7, 8 



The following lemma characterizes the coefficients qj,j-m and is useful to prove 
our next results. 



Lemma 4 Consider the coefficients Qj.j-m defined by (18). Then 

(a) qjj+i = 1, qjj < 0, ^jj-m > 0, m > 2, \mijjjj_rn = and 

Qj,j-{m+l) < Qj,j-m < Qj,j-2- 

oo 

(b) g,- = -3 + 3 X 2=^-" - 3^-". 

m=2 

oo 

m=—l 
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Proof : (a) We have that qjj+i = ajj = 1, qjj = 2^ " — 4 < 0, for 1 < a < 2 
and = 3^-° - 4 x 2^-°" + 6, which can be positive or negative depending 

on the value of a. The qjj-m, fri > 2, are of the form 

q. ._^ = (m + 2)3"" - 4(m + 1)^"° + Gm^"" - 4(m - 1)^-" + (m - 2)^-". 

Hence, 




leading to 



■ m 



3-a 



fc=4 



3 — a \ / 1 



V 



m 



- m 



-4E 

fc=4 
3-a 

1 



3 — a 
k 



m 



k CO 



k=4: 



3 — a 
k 



m 



m 



a-l 



(3 - a)(3 -a- 1)(3 - a - 2)(3 - a - 3) 24 
4! 

"(3 - a)(2 - a)(l - 24 



4! 



(35) 



Considering (35) and noting that the k odd terms of the series cancel, the 
properties (a) can be easily obtained. 

(b) In order to compute the series, let us first compute the sum of the first 
M — 1 terms. We have 



M 
m=2 



-3 + 3x2 



3— a q3— o 



3^-" + sm, 



where 



Sm = -{M - I)''"" + 3M''-" - 3(M + 1)^"° + (M + 2) 



\3-a 
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Similar to what is done in (a) we can write 



sm = M 



3-a 



3-Q 



2 / 1 \' 



■ M 



3-Q 



g 3-a. 
fc=o \ k 



2 

My 



g 3-. 

A;=0 \ A; 



~MJ 



Therefore 



fc=o \ k 



3-a- 



M 



sm = M 



3-a 



El '-"WAV 

fc=3 \ k 



Mj 



-E 

fc=3 

1 



~ / 3-a 
fc=3 \ k 



3 — a 

(3-a)(2-a)(l -a) 6 
3! JP 
■(3-a)(2-a)(l -a) 6 



3! 



M 



Clearly, we can conclude that liniM-i^oo = 0. Hence, 



M 



M 



V qj,j-m = lim V qj,j-m = -3 + 3x2 



3— a q3— a 



m=2 



m=2 



(36) 



(c) This result comes immediately from (b) and from the fact that qjj+i 
Qjj + qjj-i = 3 - 3 X 23-" + 33-". 



Remark: Note that, the previous result lead us to conclude the series, defining 
the operator (19), converges absolutely when we have a bounded function u. 



The next theorem states the method is unconditionally stable for r G [1/2, 1]. 

Theorem 5 The weighted numerical method (23) is unconditionally von Neu- 
mann stable for r G [1/2, 1] . 
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Proof: Let us insert the mode K^e*-'^ into (34). We obtain the following 



_ ^ 



a oo 



a oo 



J(j-m)e 



e'^' + (1 - r) 



The amplification factor is given by 



j—rrr~ 



K{d) 



1 - T- 



a oo 



r(4 - a) 



-imO 



m=—l 



l + (l-r) 



a oo 



r(4 - a) 



E 



-imO 



m=— 1 



Therefore |k(^)| < 1 if and only if the real part of the series is negative, that 



IS, 



^ qjj_rnCos{m9) < 0, 



m=—l 



since the imaginary part of the right side is smaller for r G [1/2, 1], because 
r > 1 — r. We can write 



E cos{m9) = {qjj+i + qjj.i) cos(^) + g^j 

m=— 1 

oo 

+ E qj,j-mCos{m9). 

m=2 

Since gjj+i + qj,j-i > 0, and qjj-m > for m > 2, 



(37) 



^ g^j.^ cos(m^) < (g^j+i + g^j.i) + g^^ + ^ g^- 

m=— 1 m=2 

Now using Lemma 3. (c), we obtain 



(3J 



^ qj,j-m cosimO) < 0. 



m=— 1 



□ 



4 Matricial form 



(39) 



We start to describe the matricial form of the numerical method, taking in 
consideration that to implement the numerical method we need to have a 
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computational bounded domain. Let us assume we consider the computational 
domain [a, b], where the mesh is defined as xj = a + jAx and we assume we 
have 

u{a,t) = 0, and u{b,t) = gh(t) given 

It is straightforward to conclude, that if u{a,t) = 0, the problem is equivalent 
to a problem defined in the whole real fine with the solution zero for x < a. 

The numerical method can be written in the matricial form 



r(4-a 



-Q U 



71+1 



r 



fx- 



r(4 - a 



-Q\W 



r(4-a) 



rb 



n+l 



1 - r)b") + p"+^ (40) 



where 
U" = 



n+T 



P 



AtTp'^+^ + (1 - r)p'l . . . AtTp^^^ 



N-l 



N-l 



contains the boundary values, /x" is a diagonal ma- 
trix with entries and Q is related to the fractional operator. The matrix 
Q = [Qj,k] has the following structure 



Q j,k 



qj,k 

Qj,j+ii k = i 
0, k>3 



l<k<i -I 
k = j 



1 
1. 



Finally the vector b" is given by 

0, 3 

assuming that Uq = and = Qbitn)- 



.N -2 



N-l. 



Remark: From Lemma 4, for qjj-i > (i.e. a > 1.5545), we can also easily 
prove our numerical method is unconditionally stable by the Gerschgorin's 
theorem applied to the iterative matrix. 



5 Numerical implementation 



The numerical experiments are carried out in two parts. First, we verify the 
accuracy and order of convergence of the numerical method to confirm the the- 
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oreticall results presented in the previous sections. Then a physical application 
is considered to reveal some of the physical phenomena, from anomalous to 
mormal diffusion. 

Consider the vectors Uapp{Ax) = {Uq, . . . , Uj\f), where Uj is the approximate 
solution, for Xj = Xq + jAx, j = 0, . . . , N at a certain time t, and Uex{A.x) = 
{u{xo,t), . . . ,u{xN,t)), where u is the exact solution. The error is defined by 
the Zoo norm as, 

\\Uex{Ax) - Uapp{Ax)\\oo = mi^^\u{xj,t) - Uj \ . (41) 



Example 1. Consider the problem with initial condition u{x, 0) = 4x^(2— a;)^, 
< X < 2 and zero otherwise. Let 

d{x) = ^r(5 - (42) 



and 

p{x, t) = -4e-*x^ [7(2 - xf + 2a{a - 7) + 6xa\ . (43) 



The exact solution is given by u{x,t) = 4e *a;^(2 — x)^, for < x < 2, and 
zero otherwise. 

In Table 1, we show the behaviour of the error (41) for different values of r 
and for At = Ax = 1/30 for the problem (42)-(43). 



T 


a = 1.2 




a = 1.4 




a = 1.5 




a = 1.8 




0.5 


4.0277x10 


-3 


3.4191x10 


-3 


3.1944x10 


-3 


2.4542x10 


-3 


0.6 


5.6194x10-^ 


-3 


4.9682x10" 


-3 


4.7877x10" 


-3 


4.2856x10" 


-3 


0.7 


7.5094x10" 


-3 


6.7573x10" 


-3 


6.5920x10" 


-3 


6.2510x10" 


-3 


0.8 


9.5429x10" 


-3 


8.6634x10" 


-3 


8.4903x10" 


'3 


8.2598x10" 


-3 


0.9 


1.1656x10" 


-2 


1.0625x10" 


-2 


1.0435x10" 


-2 


1.0283x10" 


-2 


1.0 


1.3814x10" 


-2 


1.2615x10" 


-2 


1.2403x10" 


-2 


1.2318x10" 


-2 



Table 1 

Global loo error (41) of time converged solution at t = 1 for a = 1.2, a = 1.4, 
a = 1.5, a = 1.8 and At = Ax = 1/30. 

The most accurate result is for r = 1/2. For the same problem, we observe in 
Table 2 and Table 3 that for all values of a we have second order convergence 
as expected, when r = 1/2. 
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Ax 


a = 1.2 




Rate 


a = 1.4 




Rate 


1/5 


1.5310x10" 


-1 




1.1950x10" 


-1 




1/10 


3.6239x10" 


-2 


2.0789 


3.0270x10" 


-2 


1.9811 


1/20 


9.0506x10" 


-3 


2.0015 


7.6627x10" 


-3 


1.9820 


1/40 


2.2669x10" 


-3 


1.9973 


1.9289x10" 


-3 


1.9901 



Table 2 

Global loo error (41) of time converged solution for four mesh resolutions at t = 1 
for a = 1.2, a = 1.4, At = Ax and r = 1/2. 



r 


a = 1.5 




Rate 


a = 1.8 




Rate 


1/5 


1.0884x10" 


-1 




7.9651x10" 


~2 




1/10 


2.8101x10" 


-2 


1.9535 


2.0820x10" 


-2 


1.9357 


1/20 


7.1358x10" 


-3 


1.9775 


5.4174x10" 


-3 


1.9423 


1/40 


1.8050x10" 


-3 


1.9831 


1.3974x10" 


-3 


1.9549 



Table 3 

Global loo error (41) of time converged solution for four mesh resolutions at t = 1 
for a = 1.5, a = 1.8, At = Ax and r = 1/2. 

Example 2. Consider now a second problem with initial condition u{x, 0) = 
x^, < X < 1 and boundary conditions M(0,t) = and M(l,t) = e~*. Let 

d{x) = and p{x, t) = -{l + x)e~'x\ (44) 

r(A + 1) 



The exact solution of the problem is of the form 

u{x,t) =e-^x^, xe[0,l]- (45) 



Although this problem is not defined in the whole real line we have u{0,t) = 0, 
and this can be seen as a problem for which the solution is zero when x < 0. 

In Table 4, we show the behavior of the error (41) for different weighted 
coefficients r. We observe the most accurate behaviour is again for r = 1/2. 

In Table 5 we present a comparison between our method and the methods 
presented in [16] with the same space and time steps. The second column 
shows the absolute value of the largest error calculated by the Crank- Nicolson 
scheme (before extrapolation) presented in [16] at time t = 1.0 which consists 
of assuming the fractional derivative is approximated by the shifted Griinwald- 
Letnikov formula. The third column shows the error calculated by the Crank- 
Nicolson scheme after a Richardson's extrapolation presented in [16]. The 
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T 


a = i.z 




^, 1/1 

a = i.4 




a = i.O 




a = i.O 




U.o 


b.47y2x 10 


-5 


Z.94U2X 10 


-5 


1.7o5L)X 10 


-5 


4.(J5L)9x 10 


-6 


U.D 


y.o854x lU 


-4 


7.0639X 10 


-4 


6.2104X 10 


-4 


A.bLZZX LU 


-4 


n 7 
U. i 


i.oDuy X iu 


-3 


i.ooiOX iU 


-3 


i.zzooX iU 


-3 


y.uo4o X iu 


-4 


0.8 


2.7426x10" 


'3 


2.0533x10- 


-3 


1.8233x10" 


-3 


1.3587x10" 


-3 


0.9 


3.6143x10-^ 


-3 


2.7219x10" 


-3 


2.4211x10" 


-3 


1.8110x10" 


-3 


1.0 


4.4769x10' 


-3 


3.3870x10" 


-3 


3.0166x10" 


-3 


2.2624x10" 


-3 



Table 4 

Global /oo error (41) of time converged solution for the problem (44) calculated by 
weighted numerical scheme with At = Ax = 1/30, A = 3,0 < x < 1 for different 
values of a and r. 

fourth column shows the largest absolute error for our numerical scheme with 
r = 0.5. Note that our numerical results are more accurate than the method 
given in [16]. 



Ax 


CN-GL [16] 


Extrapolated CN-GL [16] 


Weighted (r = 


: 0.5) 


1/10 


1.82265x10-3 


1.77324x10"^ 


3.5504x10" 


-5 


1/15 


1.16803x10"^ 


7.85366x10"^ 


1.6197x10" 


-5 


1/20 


8.64485x10"^ 


4.40627x10"^ 


9.1072x10" 


-6 


1/25 


6.84895x10"^ 


2.82750x10"^ 


5.8030x10" 


-6 



Table 5 

Global loo error (41) of time converged solution for the second problem calculated 
at t = 1 for the second problem with At = Ax, A = 3,0<x<l and a = 1.8. 

To conclude this example we observe the rate of convergence of the numerical 
method for different values of r 7^ 1/2. The expected convergence rate for 
r 7^ 1/2 according to section 3 is 0{At + Ax^). We consider At = Ax^ to get 
second order convergence as we observe in Table 6. 

Example 3. Finally, in order to reveal the dynamics behavior of the diffusion 
equation (1), in this example we consider equation (1) without the source 
function (which means p{x,t) = 0) on a finite domain [0,4]. We consider the 
Gaussian function 

«(x,0) = ^exp(-^^f^) 
as the initial condition, the diffusion coefficient d{x) = 1 and the boundary 
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Ay- 




(JL l.O 




xvclLt: 


r = 0.6 


1/100 


1/5 
1/10 


7.9oz5x lU 
2.1501x10" 


-4 
-4 


1.8834 




1 //I nn 
i/4UU 


1 /on 
1/ zU 


^ QVi n \/ 1 n^ 
o.o / iU X iU 


-5 


o nni 1 
z.UUii 




1/1 finn 


1 /An 

i / 4tU 


1 "^i^i 9 Y 1 n^ 

I.OOIZ A lU 


-5 


1 QQOQ 


r = 0.7 


1/25 
1/100 


1 /K 

1/5 
1/10 


1 oco^ 1 n^ 
i.zoo ( X iU 

3.5212x10" 


-3 
-4 


1.8662 




1 //inn 
i/4UU 


1 /on 
i/ZU 


Q 'yono V/ 1 n^ 


-5 


o nnoQ 
z.UUzo 




1 /I Rnn 


1 /An 

i / 4tU 


9 9n^'^ V 1 n^ 

Z.ZUOO X lU 


-5 


1 QQAR 


r = 0.8 


1/zo 
1/100 


1 /K 
1/5 

1/10 


1 'y^oQ 1 n-^ 
1. / ^^ZOX iU 

4.8915x10" 


-3 
-4 


1.8573 




1 //I nn 
i/4UU 


1 /on 

1/ ZU 


1 99nv \/ 1 n~ 
i.zzU ( X iU 


~4 


z.UU-iiD 




1 /I i^nn 

1 / IDUU 


1 /An 

i / 4tU 


ni^QA V 1 n^ 
o.uoy^t X lu 


-5 




T = 0.9 


l/zo 
1/100 


1 /c; 
1/5 

1/10 


OOP; nn \ ^ 1 n- 

z.z59Ux lU 
6.2608x10" 


-3 

-4 


1.8513 




-\ / 4 nn 
i/40U 


1 /on 
i/zO 


1 c^f^Oyi 1 n~ 


-4 


o nno/^ 
Z.UUzD 




1/1 «nn 


1 //in 


Q Q1 Q/l 1 n- 
o.y io4 X iU 


-5 


1 007"? 


T = 1.0 


1 /o 

i/zo 
1/100 


1/5 
1/10 


O ^ A QC 1 n^ 

z. ( 4oo X iU 
7.6292x10" 


-3 
-4 


1.8466 




1/400 


1/20 


1.9041x10" 


-4 


2.0024 




1/1600 


1/40 


4.7674x10" 


-5 


1.9978 



Table 6 

Global Zoo error (41) of time converged solution for the numerical scheme (23) at 
t = 1, for At = Ax^ and a = 1.8, A = 3 and different values of r. 



conditions u(0,t) = ?x(4,t) = 0. The numerical results for this example are 
calculated by the weighted scheme with r = 1/2. In this test, we take a = 0.01. 
The evolution of the non-Fickian diffusion processes for different values of a 
are given in Fig 3. The anomalous diffusion parameter exhibits the extent 
of the long tail diffusion processes of problem (1). The non-Fickian behavior 
gradually disappear when a — )■ 2. This is consistent with the experimental 
results [1,7,13,14]. Again the validity of our numerical methods is confirmed. 
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(a) a = 1.4 



(b) a = 1.6 




(c) a = 1.8 (d) a = 1.999 

Fig. 3. The evolution of u{x, t) for different anomalous diffusion coefficients a at 
different times. 



6 Conclusions 



We have derived a weighted numerical method for the fractional diffusion 
equation based on the Riemann-Liouville derivative defined in an unbounded 
domain. The numerical method is second order accurate for r = 1/2 and first 
order accurate for r G (1/2,1] because of the time discretization. We have 
proved theoretically the method converges by showing consistency and von 
Neumann stability. In the end we have presented test problems which are in 
agreement with the theoretical results. 
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